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ABSTRACT 

Using a series of three-dimensional, hydro dynamic simulations on an adaptive 
grid, we have performed a systematic study on the effect of bubble-induced motions 
on metallicity profiles in clusters of galaxies. In particular, we have studied the de- 
pendence on the bubble size and position, the recurrence times of the bubbles, the 
way these bubbles are inflated and the underlying cluster profile. We find that in hy- 
drostatic cluster models, the resulting metal distribution is very elongated along the 
direction of the bubbles. Anisotropics in the cluster or ambient motions are needed if 
the metal distribution is to be spherical. In order to parametrise the metal transport 
by bubbles, we compute effective diffusion coefficients. The diffusion coefficients in- 
ferred from our simple experiments lie at values of around ~ 10^^ cm^s~^ at a radius 
of 10 kpc. The runs modelled on the Perseus cluster yield diffusion coefficients that 
agree very well with those inferred from observations. 

Key words: 



1 INTRODUCTION 

The hot, diffuse gas that permeates clusters of galaxies, the 
intra-cluster medium (ICM), has a metallicity of about 1/3 
of the solar value. Meanwhile, the total amount of iron in the 
ICM is huge: it is larger than the total iron mass in all cluster 
galaxies. Moreover, this value does not seem to vary much 
with time, at least not up to a redshift of ~ 1 (|Tozzi et al.l 
The recent generation of X-ray observatories has pro- 
vided m aps of the radial distribution of metals in the nearby 
clusters (IPe Grandi Molendi" 2001; Fukazawa et al." 2000; 
Schmidt et al. 2002; Matsushita et al. 2002; Chu razov et al. 



20031 : iDe Grandi et al.ll2004l : iTamura et allliool ). These ob- 



servations have revealed an interesting trend: 

Galaxy clusters can be grouped into two categories de- 
pending on their X-ray surface brightness profiles: (i) clus- 
ters with a central peak in the X-ray surface brightness, 
i.e. clusters with a cool core and (ii) clusters without a 
cool core. Interestingly, these two groups show a different 
spatial distribution of metals. The clusters without cool 
cores have a nearly uniform spatial distribution of met- 
als, while clusters with cool cores have a strongly peaked 
abundance profile. Moreover, the relative abundance of el- 
ements in cool core clusters suggests that Type la super- 
novae have provided mo s t of th e metals in the central iron 
peak. iFinoguenov et al] (|2002[ ) have found that the cen- 
tral abundance peak is strongly enriched in iron while the 



bulk of the ICM outside the central region had an iron-to- 
silicon ratio close to the yields of Type II supernovae (also 
iMushotzkv et al.l ll996). This suggests that the central metal 
peak is mostly caused by Type la supernovae in the central 
cluster galaxy. However, in order to produce the large mass 
of metals in the central peak, long enrichment times are nec- 
essary. 

There are a number of observations that indicate that 
stars in massive galaxies are responsible for the metal 
enrichment of the ICM. It has been established that the 
iron mass of the ICM correlates with the optical light o f 
massive early-type galaxies in clusters ("Arnau d et al.l1l99^ 
and that the ratio of the iron mass to light and to the 
total mass is approximately constant. Thus, stellar mass 
loss from the central cluster galaxy through supernovae 
and stellar winds appears to be the prime source for the 
metals observed in the inner parts of galaxy clusters. 
However, the observed metallicity profiles are much broader 
than the stellar light profiles of the central galaxy. Hence, 
the differences in the light and metal distributions are 
interpreted as the result of transpo rt processes th a t have 
mixed the metals into the ICM. iRebusco etHI (|2005l l 
have assumed that the spreading is due to local turbulent 
motions in the ICM, and that the evolution of the metal 
density profile can be modelle d as a diffusion process. 
Using th e abundance pr o file fro m IChurazov et al. 

1 (2003) for 

Perseus, iRebusco et al.l (|2005h estimated that a diffusion 
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coefficient of the order of 2 x 10^^ cm^ is needed to 
explain the width of the observed abundance profiles. 
Based on a model of iron enrichment from stellar mass 
loss in the central clus ter galaxies of a few nearby clusters, 
iBohringer et al ] (|2004h have used the total amount of iron 
observed in the ICM to obtain constraints on the age of the 
cooling cores (i.e. the time for which the central region have 
remained relatively undisturb ed). They find very large ages 
of more than 7 Gyr (Bohrin ger et al.l ^2004) that provide 
new constraints for the modelling of the interaction regions 
and have consequences for the turbulent transport of 
metals. While it appears to be established that the metals 
produced by the central galaxy are dispersed into the ICM 
to form the broad abundance peaks, it remains unclear 
what the mechanism is via which the metals are transported. 

Currently, the most popular model that is invoked 
to explain the apparent stability of cool cores against a 
cooling catastrophe relies on heating by a central Active 
Galactic Nucleus (AGN). Radio- loud AGN are also thought 
to regulate the formation of the most massive galaxies and 
thus explain the c ut-off of the galaxy luminosity function 
at the bright end (jCroton et al.l 120061 ). Radio- loud active 
galactic nuclei drive strong outflows in the form of jets 
that inflate bubbles or lobes. The lobes are filled with 
hot plasma, an d can heat the cluster ga s in a number 
of ways fe.g. ISruggen et al.l 120 02: Briiggen Kaisei 



2002 

20041: LOmma et al 



Churazov et al, 

r-r 



Dalla Vecchia et al 



2004'; 'Basson Alexan derl l2003l: 
Revnolds et al. 2002; Ruszkowski e t al. 2004 : Briiggen et al.l 
20051 : iBrighenti Mathewsl 120061 ). 



In fact, radio sources in cooling clusters are different 
from the general population of radio sources in galaxy 
clusters. Most radio galaxies in clusters are FR I sources in 
terms of their radio power, bu t they show a great v ariety 
of morphologies and dynamics (|Owen Ledlowlll997l ). The 
central galaxy in almost every strong cooling core contains 
an active nucleus and a currently active, jet-driven radio 
galaxy. High-resolution X-ray observations of cooling flow 
clusters with Chandra have revealed a multitude of X-ray 
holes on scales < 50 kpc, often coinci dent with pat c hes o f 
radio emission. In a recent compilation, iBirzan et aD (|2004l ) 
lists 18 well documented clusters which show X-ray cavities 
with radio emission. 



~ 10% of the overall observed level of enrichment in the ICM 
within a radi us of 1.3 Mpc (s e e also iTornatore et al]l2004l ). 
Furthermore, Schindler et al.l (|2005l ) find that ram-pressure 
stripping is more efficient than quiet (i.e. non-starburst 
driven) galactic winds, at least in recent epochs since 
redshift 1. Also, they find that the expelled metals are not 
mixed immediately with the ICM, but inhomogeneities are 
visible in the metallicity maps. Relatively little theoretical 
work has gone in establishing the importance of AGN- 
induced transport of m etals in cluster centres (|Bruggenl 
I2OOI lOmma et al]l2004l ) . However, this mode is likely to be 
very important, at least at redshifts ^ 1: If AGN- induced 
motions are sufficient to quench cooling flows at the centres 
of galaxy clusters, they are similarly likely to affect the 
metal distribution in the central cluster region. Simulations 
as those presented here are essential for a proper modelling 
of the chemical evolution of clusters via semi-analytical 
techniques, as for example in iDe Lucia et al.l (|2004l ) : ICoral 
(l2C 



The broad abundance peaks in clusters have been 
produced over a time span of several gigayears. If the peaks 
have been broadened by AGN, this process has taken a 
large number of activity cycles. In order to simulate the 
transport of metals by AGN-induced flows and to capture 
the hydrodynamical details important for the interaction 
between the AGN and the ICM, a fair resolution of the 
computational mesh is necessary. This prevents us from 
simulating the bubble transport over cosmological times. 
However, we can still study the efficiency of the metal 
transport by simulating a small number of AGN cycles 
and parametrise the transport efficiency. A convenient 
parametrisa tion that has been applied previously to ob- 
servations ([Rebusco et al.l l2005l ) is based on a diffusion 
description. Even though we do not propose that the 
metal transport occurs via microscopic diffusion, transport 
processes like these may be described by a diffusion pa- 
rameter. This allows a direct comparison to observations 
without having to simulate the entire life time of the cluster. 

In this paper, we investigate the influence of buoyantly 
rising bubbles in clusters on the metal distribution in the 
ICM. In particular, we focus on the following questions: 



Buoyantly rising bubbles induce subsonic motions 
in th e cluster gas that can redistribute metals in the 
ICM (Briig genI I2OO2I ). Considering the impact of resonant 
scatterings on the surface brightness profile of the most 
prominent X-ray emission lines, evidence for gas motions 
was found in the XMM - Newton data on the Perseus cluster 
(jChurazov et al.ll2004l ). The same motions could spread the 
metals over the cluster volume although the efficiency of 
spreading could be very sensitive to the character of these 
motions. Other enrichment mechanisms include supernovae- 
driven winds, radiation- pressure- driven dust effiux (see e.g. 
lAguirre et al. _ 2001 ) and ram pressure stripping of cluste r 
galaxies ([Roediger Brueggenll2005l : iRoediger et al.lE006l ). 
Using cosmological simulations and a heur istic prescription 
for th e effect of ram-pressure stripping, IPomainko et al.l 
(|2006l ) find that ram-pressure stripping can account for 



• How far out can bubbles carry metals? 

• How efficient is transport by buoyant bubbles? 

• Does this process cause characteristic features, e.g. 
anisotropics? 

• How does the resulting metal distribution depend on 
bubble and cluster parameters? 

• What effective diffusion coefficients does this form of 
transport yield? 

• How do th ese coefficients compare with those found by 
iRebusco et aP (,2005. )? 



To this end, we simulated the effect of bubble-induced 
motions in the ICM based on the 3D hydrodynamics code 
FLASH. 
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2 METHOD 
2.1 Overview 

We model the hydro dynamical evolution of the ICM in a 
galaxy cluster. Initially, the ICM is set in hydrostatic equilib- 
rium in a static cluster potential. Details of the ICM model 
are given in Sect. 12.3.11 

We assume that a central cluster galaxy injects metals 
into the ICM with a rate proportional to its light distribu- 
tion. Details of the central galaxy are given in Sect. 12.3.21 
We now trace the distribution of the metals by injecting a 
tracer fluid into the ICM at a rate that is proportional to the 
light distribution of the central galaxy. Section [2.41 explains 
the details of the metal injection and the way the metals are 
traced throughout the simulation. 

Finally, we model the AGN activity by inflating am- 
bipolar pairs of underdense bubbles in the ICM that rise 
buoyantly and thus stir the ICM. In Section [231 we explain 
in detail how the bubbles are generated. 



2.2 Code 

The simulations we re performed with the FLASH code 
([Frvxell et al.l l2000h , a multidimensional adaptive mesh 
refinement code. FLASH is a modular block-structured 
AMR code, parallelised using the Message Passing 
Interface (MPI) library. It solves the Riemann prob- 
lem on a Cartesian grid using the Piecewise-Parabolic 
Method (PPM). The simulations are performed in 3D 
and all boundaries are reflecting. We use a simulation 

box of size (x minj 3^max) X (^/minj 2/niax) X (^^min 5 'S'max) — 

( — 145 kpc, 145 kpc)^. The cluster centre is located at 

(Xcluster,?/cluster, ^cluster) = (0,0,0). For OUr grid, WC choSC 

a block size of 16^ zones. The unrefined root grid contains 
6^ blocks. The refinement criteria are the standard density 
and pressure criteria. In the central ^ 30kpc, the grid 
is assigned a minimum refinement level of 3. Outside 
^ 100 kpc, the refinement is restricted to a maximum of 
5 levels to limit the computational effort. The resulting 
resolutions for the different runs are summarised in Table [S] 



Table 1. ICM parameters for generic runs. 



^ICM 


50 kpc 


PICMO 


10-26 gcm-3 


(3 


0.5 


TiCM 


4.7- 10'^ K 



The Piecewise-Parabolic Method uses a monotonicity con- 
straint rather than artificial viscosity to control oscillations 
near discontinuities. However, the discretisation of the 
equations still leads to a numerical viscosity. A comparison 
with the simul a tions of bubbles in a viscous medium by 
iRevnolds et al.l (|2005l ) (cf. with their figure 4) shows that 
the Reynolds number in our simulations is > 1000. 

Except by the effect of viscosity, bubbles can be sta- 
bilised by magnetic fields. It has been known for some time 
that the hot gjas in clusters o f galaxies hosts significant mag- 
netic fields ([Carilli Tavlo r 2002), with estimated typical 
magnetic field strengths of 5 x 10~^ G (5 /xG). A mag- 
netic field tangential to the fluid interface can st abilise the 
Ray leigh- Taylor instability (|Chandrasekharlll96lh . provided 
the restoring tension generated by bending of the field lines 
exceeds the buoyancy force driving the instability. The sta- 
bilising effect of intracluster magnetic fields was pointed 
out by De Young (2003), who derived analytic conditions 
for the stabilisation of relic radio bubbles in the hot ICM. 
These calculations showed that the field strengths observed 
in many clusters will stabilise the bubble interface in the 
Unear regime, and hence the cluster magnetic fields could 
account for why the relic radio bubbles are seen as in- 
tact objects at such late times. Additional support for this 
idea i s found in th e two-dimensional MHD calculations of 
Briigg ^n" Kaiserl (^2001), who applied a /3 ^ 10 field inside 
the bubbles that was aligned with the bubble surface. While 
viscosity and magnetic fields can affect the dynamics of the 
bubbles and thus also the resultant transport of metals, the 
effects are still poorly understood. Here we only study the 
case of inviscid and unmagnetised fluids. 



There have been some sugges tions that the ICM 
may have a non- negligible viscosity (Reynolds et al. 1 l2005l : 
iRuszkowski et al . 2004). However, the effective viscosity of 
the ICM remains unknown because the physics of such dilute 
and magnetised plasmas is poorly constrained. In particular, 
the effect of magnetic fields on the macroscopic viscosity is 
unclear. Even minute magnetic fields lead to small proton 
gyroradii and are thus expected to suppress the viscosity 
efficiently. However, it has been pointed out that the ex- 
ponential divergence of neigbouring field lines in a tangled 
magnetic field may lead to only a very modest suppression 
of the viscosity. iRevnolds et al.l (|2005r ) have simulated the 
buoyant rise of bubbles in a medium with a kinematic vis- 
cosity. They found that even a modest shear viscosity (cor- 
responding to 1/4 of the Spitzer value) can quench fluid in- 
stabilities and keep the bubbles intact (especially see figure 
4 in their paper). Other simulati ons of buoyant bubbles in 
a viscous ICM were performed by ISiiacki &: Springel (|2006l l 
using smoothed particle hydrodynamics. 

FLASH solves the equations of inviscid hydrodynamics. 



2.3 Model cluster 

2. 3. 1 ICM distribution 

In this paper, we study two different cluster models. In the 
first, generic one, we take a constant ICM temperature and 
the pressure is set by assuming hydrostatic equilibrium. Ta- 
ble [T] lists the ICM parameters. 

The density of the ICM follows a /3-profile, 



PlCM(r) — piCMO 



1 + 



-3/2/3 



(1) 



In addition to the generic cluster model, we set up 
a cluster modelled on the brightest X-ray cluster A426 
(Perseus) that has been studied extensively with Chandra 
and XMM-Newton. 

The electron density Ue and the temperature Te profiles 
used here are based on the deprojected XMM-Newton data 
(Churazov et al. 2003, 2004) which are also in broad agree- 
ment with the ASCA (Allen & Fabian 1998), Beppo-Sax (De 
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Table 2. Central galaxy - metal injection parameters. 



a 10 kpc 

M^etaio 952M0Myr-i 



Grandi & Molendi 2001, 2002) and Chandra (Schmidt et al. 
2002, Sanders et al. 2004) data. Namely: 



metal injection is arbitrary, fp and / can be scaled to the 
total metal injection rate of the central galaxy, Mmetaio- In 
the numerical implementation, the metal fraction at each 
position X is updated according to 



in each timestep. 



fold{x)p{x) + y6nietal(x) • At 

p{x) 



(7) 



4.6 X 10"^ 



4.8 X 10" 



[1 + (J_)2]1.8+ [l + (_g_)2]0, 

and 



87 



cm 



Te = 7 X 



[1 + (7t)^1 
[2.3 + (^)3 



keV, 



(2) 



(3) 



where r is measured in kpc. The hydrogen number den- 
sity is assumed to be related to the electron number density 
as riH — ne/1.2. 

2.3.2 Central galaxy 

The metal injection rate in the central galaxy is assumed to 
be proportional to the light distribution. It is modelled with 
a Hernquist profile given by: 



Pmetal(r) — 



MmetalO « 
271 



(r + a)3 



(4) 



Note that the injection rate could also be time-dependent. 
iRebusco et al.l ([2005) have used a time-dependent metal in- 
jection rate that accounts for the higher supernova rate in 
th e past and the evo lution of the stellar population (also 
se enzini et al.lll993h . However, we follow the evolution of 
th e cluster for on l y abo ut 1 Gyr. For typical cases studied 
in iRebusco et all (|2005 l). the metal injection rate does not 
change significantly over this time. Hence in our simulations 
we use a temporally constant metal injection rate. Table [2] 
lists the parameters for the metal injection. 



2.4 Tracing the metals 

The FLASH code offers the opportunity to advect mass 
scalars along with the gas density. In order to be able to trace 
the metal distribution, we utilise one of the mass scalars, /, 
to represent the local metal fraction in each cell: 

Pmetal 



/ 



(5) 



Hence, the quantity fp gives the local metal density, pmetai- 
The metal density obeys the continuity equation including 
the metal source: 



d 



(6) 



In order to model the metal injection by the central galaxy, 
we assume that the metal fraction is small at all times. 
Hence, we can neglect pmetai as a source term in the conti- 
nuity equation for the gas density, p. As a consequence, the 
system of hydrodynamical equations including Eq.[6]is iden- 
tical for the set of variables (pmetai, /) and [A • pmetai, ^ • /), 
where A is a constant. In other words, the amplitude of the 



2.5 Bubble generation 

Bubbles in the ICM are thought to be inflated by a pair of 
ambipolar jets from an AGN in the central galaxy. In this 
simple picture, the jets inject energy into small volumina 
of the ICM at their terminal points. Thus, this region is 
in overpressure and expands until it reaches pressure equi- 
librium with the surrounding ICM. The result is a pair of 
underdense, hot bubbles. 

In our simulations, we generate pairs of bubbles on op- 
posite sides of the cluster centre. There are various ways in 
which one could produce these bubbles and we have tested 
the following two methods: 

• In the first method, we infiate the bubbles by injecting 
energy into two small spheres of radius rinj with a constant 
rate e (in ergcm~^s) over an interval of length Tinj. The 
total injected energy is ^inj = e47rrfnjrinj per bubble. The 
gas inside these spheres is heated and expands similar to 
a Sedov explosion in a few Myr to form a pair of bubbles. 
The parameters rinj, ^inj and Tinj are chosen such that the 
resulting bubbles show approximately a density contrast of 
ph/pamh = 0.03 to the surrouuding ICM and approximately 
a radius of rbbi- The expansion is much shorter than the rise 
of the generated bubbles. In addition to the bubbles, the 
explosion sets off shock waves that move through the ICM. 
Incidentally, the dependence of the bubble dynamics on the 
density contrast, pb/pamb, is weak provided that pb/pamb ^ 
1. 

• In the second method, we evacuate the bubble regions 
by removing gas while keeping them in pressure equilibrium 
with their surroundings. This evacuation method does not 
produce the shock waves generated by the infiation method. 
We implement gas mass sinks inside two spheres of radius 
rbbi- Inside these regions, gas is removed with a certain rate 
p over a time Tevac- The evacuation is not done instanta- 
neously to prevent numerical problems. However, 

'T'evac IS 

small compared to the evolution timescale of the bubbles. 
The sink rate p is set to decrease the the density inside 
the bubbles down to a density contrast of pb/pamb com- 
pared to the surrounding ICM. In order to conserve the 
metal mass, the metal fraction also needs to be updated dur- 
ing the evacuation. The evacuation method does not cause 
shocks that move through the ICM and is, hence, computa- 
tionally cheaper. The evacuation method models the phase 
of the evolution after the radio-loud AGN has infiated a 
low- density bubble which has exp anded to achieve pressur e 
equilibrium with its surroundings (Briiggen Kaiserl l200lh . 
The infiation of the lobes itself cannot be modelled entirely 
self-consistently as one would have to resolve scales all the 
way down to the accretion disk of the AGN. However, jets 
are observed to infiate spherical bubbles that are in near 
pressure equilibrium with the ICM. In order to model the 
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Figure 1. Evolution of the metal density profile in the case with- 
out bubbles. The metal density increases with time proportional 
to the light profile of the central galaxy. The injection rate is 
assumed to be constant in time. It is apparent that numerical 
mixing is very small. The noise increases at larger radii as the 
refinement levels go down. 

buoyant rise of these AGN-blown bubbles, it is customary 
to start with an underdense, spherical bubble. Since gas is 
removed from a small volume in the process of forming the 
bubble, mass is not conserved in this method. However, this 
does not affect the density profile of the cluster on the time 
scales that we consider here. In Fig. [8] we show that the dif- 
ferences between methods EVAC and INFL are very small. 

For both methods, we generate a pair of bubbles every 
Tbbi- For every bubble generation cycle, the position of the 
generated bubbles changes. For the next pair of bubbles, the 
axis that connected the previous bubbles is rotated around 
the 2/- axis. The angle of rotation is fixed by the demand that 
after rirotbbi cycles the bubbles are again inflated in the first 
position. The parameters for our runs are summarised in 
Table El 



3.2 Single pair of bubbles 

In the next step, we produced a single pair of bubbles by 
evacuating two spherical regions of radius 10 kpc off-centre 
and on opposite sides of the cluster centre. The distance 
of the bubble centre from the cluster centre is 16 kpc. We 
studied the effect of this single pair of bubbles on the metal 
profile in the cluster. In particular, we wanted to compute 
how much mixing a single pair of bubbles can cause and 
what the long term fate of these bubbles is. Figures [2] to [3] 
show the results from this run. 

As found in previous studies, the bubbles start to 
rise buoyantly and fragment. While rising, they uplift the 
metals from the centre and distribute them in their wake. 
Quite quickly the bubble fragments into a torus. Associated 
with this torus is a vortical flow that causes gas to ascend 
through the axis of the torus. Along the axis of the torus, 
metals seem to be transported rather fast. However, this 
metal spreading seems to be restricted to a narrow region 
around this direction. Spreading along other directions 
seems negligible. The effect is that the metals are spread 
out significantly along the axis of the bubbles leading to a 
very elongated metal distribution. Globally, the effect of a 
single pair of bubbles is not very significant (see cumulative 
metal mass profile in Fig. [3]). 

The bubbles do not continue to rise indefinitely and 
even after 1 Gyr they do not rise beyond 100 kpc from the 
centre. When they fragment, their rise velocity decreases. 
This can easily be seen by equating the drag force to the 
buoyance force, i.e. 

-7rr pg^ -Cr pv , (8) 

where v is the terminal velocity, p the density of the ambient 
medium, g the gravitational acceleration, r the radius of 
the bubble and C a constant of order unity. This yields the 
terminal velocity 



3 RESULTS 

Table Ogives an overview of the simulations described in this 
paper. Primarily, we test the dependence of the resulting 
metal distribution on the bubble position, size and recur- 
rence time, the background model, the method with which 
bubbles are produced and the resolution of the computa- 
tional grid. 

3.1 Stability test 

Before stirring the ICM with an AGN, we tested whether 
the cluster remains in hydrostatic equilibrium for a typical 
simulation time of about 1 Gyr. To this end, we ran a sim- 
ulation where we injected metals but did not produce any 
bubbles. We found that, both, gas density and metal density 
profiles remained stable for at least 1 Gyr. Figure [1] shows 
the evolution of the metal density profile. 

Even after 1 Gyr the 3 refinement levels can be seen in 
the profiles (where the "line" becomes thicker). Numerical 
velocities are of the order of 1 kms~^ inside the grid and of 
15kms~^ at the grid boundary. 



Thus, the rise velocity oc ^/rbbi (see e.g. lChurazov et aP 

(|200lh V The fragments form a mushroom-shaped plume, not 
unlike the radio structures observed in M87, and this struc- 
ture covers a larger azimutal angle. There could be a certain 
radius where most metals are deposited. 

Two effects can lead to a more isotropic mixing of met- 
als: (i) ambient motions in the cluster that are triggered 
by galaxy motions and merger activity and (ii) a series of 
bubbles that are launched at changing locations in the clus- 
ter, for example by a precessing jet or by inhomogeneities 
of the accretion flow near the jet. In the next section, we 
investigate the second scenario. 



3.3 Series of bubbles 

In the next set of simulations, we look at the effect of a suc- 
cession of bubbles that are generated at changing positions 
inside the cluster. 
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Table 3. Simulation runs. "•" means same value as in previous line. "/" means value irrelevant for this run. The first column gives 
the coordinates of the initial bubbles centres, rbbl is the radius of the bubble, Tbbi is the recurrence time of the bubbles, E'inj the total 
injected energy, rirotbbl the number of precessions of the bubble, {x / y / Zmin / m&^) the size of the computational domain, and Ax the side 
length of the smallest computational cell. 



name 


(a^bbbybbb ^bbi) 
in (kpc, kpc, kpc) 


^bbl 
/kpc 


'''bbl 


^inj 

/Myr 


^rotbbl 

/erg 


{x / y / Ztcciiyi I max) 


effective 
no of cells 


Ax 
(kpc) 


SINGLE_BBL 


(10,12,0) 


10 


/ 


/ 


/ 


±146 


768^ 


0.38 


EVAC 
EVAC_HR 






50 


/ 

/ 


4 




192^ 
384^ 


1.52 
0.76 


INFL 
INFL_HR 








6.68 • 10^'^ 






192^ 
384^ 


1.52 
0.76 


TAUBBL200 






200 


evac 




±340 


448^ 


1.52 


DISTCTR 
LARGEBBL 


(20,24,0) 


20 


50 













PERSEUS_50 (10,12,0) 10 50 / 4 ±146 192^ 1.52 

PERSEUS_200 (10,12,0) 10 200 / 4 ±146 192^ 1.52 



S.S.I Distance of bubbles from cluster centre 

One would expect that the starting position of the bubbles 
with respect to the central galaxy affects the amount of met- 
als that can be upUfted from the galaxy. This is indeed the 
case as we see by comparing runs EVAC and DISTCTR. If 
bubbles start too far out, they cannot uplift a substantial 
amount of metals because most metals are injected mainly 
in the cluster centre (steep metal injection profile) as can be 
seen from Fig. 2] 

In order to compare metal fraction, metal density and 
cumulative metal mass at different times, these quantities 
are normalised to Mmetaio • t. From this we can conclude 
that the initial position of the bubbles with respect to the 
scale radius of the metal injection rate has an impact on the 
mixing efficiency of the bubbles. However, as the bubbles 
are produced by jets that are launched from supermassive 
black holes at the centre of the galaxy, the bubbles will never 
be inflated too far from the stars that produce the metals. 
While in principle the mixing efficiency depends on the ini- 
tial position of the bubble, the nature of the bubbles will 
ensure that this will not cause great differences in the metal 
profiles. Our generic run has a distance of the bubble centre 
from the core of the cluster of ^ 16 kpc. In reality, bubbles 
occur at a broad range of distances from the centre. The 
list of known bubbles in iDunn et al.l (|2005h shows that the 
distances that we probe lie well within the observed ones. 



S.S.2 Bubble size 

Obviously, one important parameter is the size of the bub- 
bles. In our standard scenario we chose a bubble radius of 10 
kpc which is motivated by the sizes of some of the cavities 
observed in the X-ray surface brightness maps. For example, 
the cavities of the inner bubbles in Perseus have a radius of 
13 kpc. To test the dependence of the transport efficiency on 
the bubble size, we compared our standard run with a bub- 
ble of twice its radius (8 times its volume). See Fig.[5]for the 
cumulative metal mass profiles. The main inference is that 
larger bubbles can carry metals out to larger distances. As 



can be seen in Fig. O after about 500 Myrs the large bubble 
(r = 20 kpc) carries the metals out to a distance of about 
50 kpc, while the 10 kpc bubble only lifts them to about 25 
kpc. 

The same trend is apparent from Fig. [6] which shows 
the normalised radial metal density profile for three differ- 
ent times. The main effect of the bubbles can be seen in 
the points elevated above the original line. In these cells, 
the metal density is enhanced because the bubbles have car- 
ried the metals to larger radii. It is evident that the larger 
bubbles are more efficient at transporting metals outwards. 

This behaviour is not surprising because larger bubbles 
rise faster as we saw in Eq. (|9]). As the bubbles rise faster, 
they accelerate the metals in its wake to higher velocities 
and they travel further out until instabilities have broken 
them up such that they do not travel much further. 



S. S. S Bubble frequency 

Clearly, the efficiency with which bubbles can transport met- 
als depends on the frequency with which bubbles are pro- 
duced by the AGN. In the previous simulations, we employed 
a recurrence time of 50 Myr. For comparison, we ran a sim- 
ulation with a recurrence time of 200 Myrs, all other param- 
eters being the same as in run EVAC. 

The difference in, both, radius up to which metals are 
carried and degree of deformation of cumulative metal mass 
profile can be seen in Fig. [71 After 800 Myrs, comparing the 
blue and grey lines reveals the difference between these two 
runs. The runs with a longer recurrence time leads to ap- 
preciably less mixing. A recurrence time that is four times 
longer, roughly leads to only about 1/3 as much mixing 
by the end of our simulation. Thus, the degree of mixing 
scales roughly with the frequency of bubble generation as 
one would expect. 
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Figure 2. Simulation of a single bubble pair: Slice through the centre of the computational box showing the local gas density (left 
column) and the metal fraction (right column) at 60 Myr (top) and 250 Myr (bottom). The grid lines indicate the boundaries of the 
blocks. Each block contains 8^ computational cells. 



3.3.4 Comparison between evacuation and inflation 
method 

In Fig. [HI we investigate the sensitivity of the metal profiles 
to the way we produce the bubbles. We generate the bubbles 
by either injecting hot gas into a number of computational 
cells or by removing gas mass from a spherical region of 
space while keeping the pressure constant. Reassuringly, the 
differences between the two methods are small. In the evac- 
uation method the bubbles seem to be able to carry metals 



out a little bit further. Also the global effect (flattening of 
cumulative metal mass profile, see Fig. [8} is a bit stronger 
in the evacuation scenario. However, all properties concern- 
ing the transport of metals in the cluster are very similar 
in both runs. The evacuation method is computationally 
slightly cheaper because the timesteps during the inflation 
phase become quite small. Projections of the gas density, the 
X-ray surface brightness and the normalised metal density 
for different timesteps are shown in Fig. [9] Therefore, we 
produce the bubble by evacuation in all our other runs. 
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Figure 3. Normalised cumulative metal mass for the run with a single pair of bubbles. These masses are averages in shells centred 
on the cluster. The thick solid line shows the cumulative metal mass without mixing. The bottom panel shows the difference between 
cumulative mass with and without mixing. 



3.3.5 Resolution test 

In these kinds of simulations, it is essential to assess the er- 
rors caused by the spatial discretisation. We can compare 
runs EVAC and EVAC_HR which have an effective resolu- 
tion of 192^ and 384^ zones respectively. The morphologies 
of the bubbles are very similar up to 250 Myr and the cu- 
mulative mass profiles are practically indistinguishable. In 
Fig. [TOl we show the normalised cumulative metal mass pro- 
files for runs with resolutions of 192^ zones and 384^ zones. 
The resultant metal profiles suggest that the simulation have 
converged numerically. 



3.4 The case of Perseus 

Having discussed a somewhat generic cluster model above, 
we now turn to the brightest X-ray cluster A426 (Perseus) 
that has been studied extensively with Chandra and 
X MM- Newton. Our simulations of the Perseus cluster show 
two things: (i) In the case of Perseus, the metal distribution 



is somewhat broader. This is a result of the greater pres- 
sure gradient in Perseus which leads to a greater buoyancy 
force and thus a higher velocity of the bubbles, (ii) We find 
that despite varying the initial positions of the bubbles, the 
resulting metal distribution remains fairly elongated and de- 
cidedly non-spherical. 

Fig. II II shows the local metal fraction at three different 
times in slices through the centre of the cluster. The left 
column is for the generic background cluster (EVAC) and 
right column for the simulation based on Perseus. For the 
Perseus model, the metal disribution is somewhat broader 
than in our generic model. However, all qualitative features 
discussed above for the generic model remain valid for the 
cluster based on Perseus. 



3.5 Estimate of diffusion constant 

In this section, we wish to parametrise the transport of met- 
als by computing an effective diffusion constant. We will 
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Figure 4. Evolution of the normalised metal density profile (top panel) and the normalised cumulative metal mass profile (metal mass 
inside r, middle panel). The bottom panel shows the difference between mixing and no mixing. This figure shows the impact of the 
launch position of the bubbles. We compare runs EVAC and DISTCTR. 



begin by defining a cumulative metal mass profile 

MmetAR)^ J PmetA^dV (10) 

r<R 

and the averaged metal density profile 

Pmetal(r) = — . (11) 

On the basis of these quantities, we c an estimate a d i ffusio n 
constant that can be compared with iRebusco etHI (|2005l ). 
We start from the spherically symmetrical diffusion equation 
including the metal source term 

^Pmetal (r,t) 



dt 



d d _ 

— D — Pnietal(r, t) + /C>metal (r) , (12) 



where D is the diffusion constant, and Pmetai(^) is the metal 
source profile according to Eq. 21 Discretising this between 
the radii Ri and R2 (i.e. AR — R2 — Ri and R = {Ri + 
i^2)/2) and timesteps ti and t2 (At — t2 — ti and t — {ti -\- 
12) 12) yields 



Pmetal Pmetal 



At 



D{R,t) 



^Pmetal(i?2, t) 



— pmetal{R) + 
Pmetal 



AR 



(13) 



where we have assumed that D varies only slowly with ra- 
dius, i.e. D{Ri,t) ^ D{R2,t) ^ D{R,t). Reordering yields 



D{R,t) = 

Pmetal {R,t2) Pmetal (i?,tl) 



AR 



- Pmetal(i^) At 



At [^Pnaetal(i?2,t) - |:Pnaetal(i?l,t)] 



(14) 



For the runs described in the previous sections, we have 
computed the radial profiles of the diffusion constants in 

Fig.im 



4 DISCUSSION 

In our simulations, the dynamics of the ICM is entirely 
driven by underdense bubbles that are initially at rest and 
then start to rise by the action of buoyancy. This is clearly 
a simplified picture. In reality, these bubble are inflated by 
relativistic jets whose physics is still very poorly understood. 
Real AGN jets have a dynamics that can only be fully cap- 
tured with special-relativistic MHD codes that need to re- 
solve scales much below the scales resolved in our simula- 
tions. 

However, bubbles such as those set up in our sim- 
ulations are similar to the bubbles observed in galaxy 
clusters such as the Perseus cluster. The bubbles in Perseus 
and other clusters are nearly in pressure equilibrium with 
the ambient medium, they are nearly spherical and move 
subsonically through the ICM. Thus, we ignore the entire 
process that has led to the formation of these bubbles and 
focus solely on the buoyancy- driven motions of the bubbles. 

There have been some attempts to simulate the 
inflation of the bubbles by a j et an d assess its consequences 
for the ICM. lOmma et all (|2004l l simulate the effect of 
slo w jets on a cool core in a galaxy cluster. Simulations 
by IVernaleo &: Revnoldsl (|2006h show that jets launched in 
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Figure 5. Evolution of the normalised metal density profile (top panel) and the normalised cumulative metal mass profile (metal mass 
inside r, middle panel). The bottom panel shows the difference between mixing and no mixing. Here we illustrate the impact of the initial 
bubble position and the initial bubble sizes. Shown are results from runs EVAC, DISTCTR and LARGEBBL. 



a hydrostatic cluster fail to heat the cluster on long time 
scales because the jet forms a low-density channel through 
which the material can flow freely, carrying its energy out 
of the cooling core. 

If bubbles in a hydrostatic cluster are launched repeat- 
edly at the same position, the resulting metal distribution 
becomes very aspherical and elongated along the axis of the 
bubble. If, however the bubbles are launched at changing 
positions, the metal distribution is somewhat isotropised. 
However, as is apparent e.g. from Fig. [11] even when we vary 
the bubble position within a cone of ^ 30deg, the result- 
ing metal distribution remains elongated and non-spherical. 
Thus, if the metal distribution should turn out to be spher- 
ical, this would necessitate additional motions caused for 
example by mergers or by stirring by galaxies in order to 
isotropise the metals. 

In real clusters, the ICM is not in hydrostatic equi- 
librium and can show strong motions. Galaxies that fly 
through the cluster produce motions and turbulence. 
Merger activity can stir the IGM even more violently. 
These motions will all contribute to dispersing the metals 
throughout the ICM. Moreover, density anisotropics can 
cause the bubbles to move away from the axis of the AGN 
jet, as the bubbles follow the local density gradient. Indeed, 
in quite a few clusters, most famously in M87, the bubbles 
are widely dispersed through the cluster and d o not trace 
the alleged axis of the AGN jet. iHeinz et al.1 (|2006h have 
performed three-dimensional adaptive-mesh simulations of 



jets situated at the centre of a cluster drawn from a cos- 
mological simulation. These simulations show that cluster 
inhomogeneities and large-scale flows have a significant 
impact on the morphology of the bubbles. Moreover, these 
motions help to prevent the format ion of channe l s alon g 
the bubbles axes as found bv .Vernaleo Revnoldsl (|2006l ). 

The diffusion coefficients inferred from our simple ex- 
periments lie at values of around ^ 10^^ cm^s~^ at a radius 
of 10 kpc. Right at the centre, they rise to about ^ 10^^ 
cm^s~^ and fall to about ^ 10^^ cm^s~^ at a radius of 100 
kpc. Modulo a factor of a few, the coefficients for all runs lie 
around these values. Differences are induced by the bubble 
sizes, the initial positions, the recurrence times and the pres- 
sure profile in the cluster. Interestingly, the runs modelled 
on the Perseus cluster yield diffusion coefficients that agree 
very well with t hose estimated by ass uming a simple diffu- 
sion model fsee iRebusco et aD (|2005l \V As we do not take 
into account motions that may have been induced by merg- 
ers etc., the values for the diffusion coefficients constitute 
lower limits on the real values. Thus, we can conclude that 
AGN-induced motions are sufficient to explain the broad 
abundance peaks in clusters. It is also striking that the dif- 
fusion coefficients fall so steeply with radial distance from 
the centre. Over two decades in radius they decrease by al- 
most 10 orders of magnitude. Very approximately, they fall 
like . This strong radial dependence of the transport effi- 
ciency is, partially, a result of the three-dimensional nature 
of the bubble- induced motions, i.e. uplifted metals are di- 
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Figure 6. Normalised radial metal density profile for three different times. A random subset of cells is chosen (same weight for each 
radius). Metal density is normalised to tMmetalO- Left column is for EVAC, right column for LARGEBBL. The main effect of the bubbles 
can be seen in the points elevated above the original line. In these cells, the metal density is enhanced because the bubbles have carried 
the metals to larger radii. It is apparent that the larger bubbles are more efficient at transporting metals outwards. 



luted over the entire radial shell. A second effect is the de- 
caying "lift" power of the bubbles as they rise. The bubbles 
accelerate only at the beginning of their ascent through the 
cluster. However, very quickly instabilities set in that slow 
down the bubble until it finally comes to a halt at about 
100 kpc from the centre. The maximal lift power is attained 
when the bubbles are still intact and have their largest veloc- 
ity. This occurs at radii of about 20 kpc - 30 kpc, and this co- 
incides with peaks in the diffusion coefficient as can be seen 
from Fig. [121 particularly for runs EVAC and PERSEUS_50. 

Finally, we may add that if hydrodynamic instabilities 
are suppressed by magnetic fields or viscous forces, for ex- 
ample, the bubbles would stay intact for longer and could 
rise to greater distances thus also enhancing the diffusion 
coefficients at larger radii. 



5 SUMMARY 

In a series of numerical experiments, we have performed a 
systematic study on the effect of bubble-induced motions 
on metallicity profiles in clusters of galaxies. In particular, 
we have studied the dependence on the bubble size and po- 
sitions, the recurrence times of the bubbles, the way these 
bubbles are inflated and the underlying cluster profile. 

On the basis of our 3D hydro simulations, we wish to 
make the following main points: 

• Larger bubbles lead to mixing to larger radii. For re- 
alistic parameters and in the absence of stabilising forces, 
bubbles hardly move beyond distances of 150 kpc from the 
centre. 

• Mixing scales roughly with bubble frequency. 
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Figure 7. Evolution of normalised metal density profile (top panel) and and normalised cumulative metal mass profile (metal mass 
inside r, middle panel). The bottom panel shows the difference between mixing and no mixing. Comparison of two runs with different 
bubble evacuation timescale (r^bi) (runs EVAC and TAUBBL200). Different line styles and colours correspond to different times and 
runs, see legend. The black line is the prediction of what should happen without any mixing. 



• The metal distribution is not sensitive to the way of 
how bubbles are inflated. 

• In hydrostatic cluster models, the resulting metal dis- 
tribution is very elongated along the direction of the bub- 
bles. Anisotropics in the cluster and/or ambient motions are 
needed if the metal distribution is to be spherical. 

• The diffusion coefficients inferred from our simple ex- 
periments lie at values of around ^ 10^^ cm^s~^ at a radius 
of 10 kpc. 

• The runs modelled on the Perseus cluster yield diffusion 
coefficients that a^ee very wel l with those inferred from 
observations (see iRebusco etHI (|2005l )). 
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Figure 9. Projections of gas density (left column), (middle column) and normalised metal density (right column) for different 
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Figure 10. Left: Normalised cumulative metal mass profile (metal mass inside r) for run EVAC with an effective resolution of 192^ 
zones (top) and for run EVAC_HR with an effective resolution of 384^ zones (bottom). The upper line in both plots corresponds to the 
case without mixing. Right: Radial metal density profile normalised to tMmetalO for the 192^ run (top) and the 384^ run (bottom). All 
data corresponds to a time of 400 Myrs after the start of the AGN activity. 
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Figure 11. Slice through the computational box showing local metal fraction at three different times. Left column is for the generic 
background cluster (EVAC), and right column for the simulation based on Perseus (PERSEUS_50). 
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Figure 12. Diffusion coefficients for the runs EVAC (top left), INFL (top right), TAUBBL200, (middle left), DISTCTR (middle right), 
PERSEUS_50 (bottom left) and PERSEUS_200 (bottom right). The different curves denote different times. 



